home *** CD-ROM | disk | FTP | other *** search
/ IRIX 5.3 for Indy R4400 / IRIX 5.3 for Indy R4400 175MHz.img / dist / eoe2.idb / usr / src / rcs / DCsource.Z / DCsource / decimal.c < prev    next >
C/C++ Source or Header  |  1995-02-28  |  34KB  |  1,407 lines
  1. /* 
  2.  * Arbitrary precision decimal arithmetic.
  3.  *
  4.  * Copyright (C) 1984 Free Software Foundation, Inc.
  5.  *
  6.  * This program is free software; you can redistribute it and/or modify
  7.  * it under the terms of the GNU General Public License as published by
  8.  * the Free Software Foundation; either version 2, or (at your option)
  9.  * any later version.
  10.  *
  11.  * This program is distributed in the hope that it will be useful,
  12.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14.  * GNU General Public License for more details.
  15.  *
  16.  * You should have received a copy of the GNU General Public License
  17.  * along with this program; if not, you can either send email to this
  18.  * program's author (see below) or write to: The Free Software Foundation,
  19.  * Inc.; 675 Mass Ave. Cambridge, MA 02139, USA.
  20.  */
  21.  
  22. /* Some known problems:
  23.  
  24.     There are some problems in the behavior of the decimal package
  25.     related to printing and parsing.  The
  26.     printer is weird about very large output radices, tending to want
  27.     to output single ASCII characters for any and all digits (even
  28.     in radices > 127).  The UNIX bc approach is to print digit groups
  29.     separated by spaces.  There is a rather overwrought workaround in
  30.     the function decputc() in bcmisc.c, but it would be better if
  31.     decimal.c got a fix for this.  */
  32.  
  33. /* For stand-alone testing, compile with -DTEST.
  34.    This DTESTable feature defines a `main' function
  35.    which is a simple loop that accepts input of the form
  36.    number space op space number newline
  37.    where op is +, -, *, /, %, p or r,
  38.    and performs the operation and prints the operands and result.
  39.    `p' means print the first number in the radix spec'd by the second.
  40.    `r' means read the first one in the radix specified by the second
  41.    (and print the result in decimal).
  42.    Divide in this test keeps three fraction digits. */
  43.  
  44. #include "decimal.h"
  45.  
  46. #define MAX(a, b) (((a) > (b) ? (a) : (b)))
  47. #define MIN(a, b) (((a) > (b) ? (b) : (a)))
  48.  
  49. /* Some constant decimal numbers */
  50.  
  51. struct decimal decimal_zero = {0, 0, 0, 0, 0};
  52.  
  53. struct decimal decimal_one = {0, 0, 1, 0, 1};
  54.  
  55. /*** Assumes RADIX is even ***/
  56. struct decimal decimal_half = {0, 1, 0, 0, RADIX / 2};
  57.  
  58. decimal static decimal_add1 (), decimal_sub1 ();
  59. static void add_scaled ();
  60. static int subtract_scaled ();
  61.  
  62. /* Create and return a decimal number that has `before' digits before
  63.    the decimal point and `after' digits after.  The digits themselves are
  64.    initialized to zero.  */
  65.  
  66. decimal
  67. make_decimal (before, after)
  68.      int before, after;
  69. {
  70.   decimal result;
  71.   if (before >= 1<<16)
  72.     {
  73.       decimal_error ("%d too many decimal digits", before);
  74.       return 0;
  75.     }
  76.   if (after >= 1<<15)
  77.     {
  78.       decimal_error ("%d too many decimal digits", after);
  79.       return 0;
  80.     }
  81.   result = (decimal) malloc (sizeof (struct decimal) + before + after - 1);
  82.   result->sign = 0;
  83.   result->before = before;
  84.   result->after = after;
  85.   result->refcnt = 0;
  86.   bzero (result->contents, before + after);
  87.   return result;
  88. }
  89.  
  90. /* Create a copy of the decimal number `b' and return it.  */
  91.  
  92. decimal
  93. decimal_copy (b)
  94.      decimal b;
  95. {
  96.   decimal result = make_decimal (b->before, b->after);
  97.   bcopy (b->contents, result->contents, LENGTH(b));
  98.   result->sign = b->sign;
  99.   return result;
  100. }
  101.  
  102. /* Copy a decimal number `b' but extend or truncate to exactly
  103.    `digits' fraction digits. */
  104.  
  105. static decimal
  106. decimal_copy_1 (b, digits)
  107.      decimal b;
  108.      int digits;
  109. {
  110.   if (digits > b->after)
  111.     {
  112.       decimal result = make_decimal (b->before, digits);
  113.       bcopy (b->contents, result->contents + (digits - (int) b->after), LENGTH(b));
  114.       return result;
  115.     }
  116.   else
  117.     return decimal_round_digits (b, digits);
  118. }
  119.  
  120. /* flush specified number `digits' of trailing fraction digits,
  121.    and flush any trailing fraction zero digits exposed after they are gone.
  122.    The number `b' is actually modified; no new storage is allocated.
  123.    That is why this is not global.  */
  124.  
  125. static void
  126. flush_trailing_digits (b, digits)
  127.      decimal b;
  128.      int digits;
  129. {
  130.   int flush = digits;
  131.   int maxdig = b->after;
  132.  
  133.   while (flush < maxdig && !b->contents [flush])
  134.     flush++;
  135.  
  136.   if (flush)
  137.     {
  138.       int i;
  139.  
  140.       b->after -= flush;
  141.       for (i = 0; i < LENGTH (b); i++)
  142.     b->contents[i] = b->contents[flush + i];
  143.     }
  144.  
  145. }
  146.  
  147. /* Return nonzero integer if the value of decimal number `b' is zero.  */
  148.  
  149. int
  150. decimal_zerop (b)
  151.      decimal b;
  152. {
  153.   return !LENGTH(b);
  154. }
  155.  
  156. /* Compare two decimal numbers arithmetically.
  157.    The value is < 0 if b1 < b2, > 0 if b1 > b2, 0 if b1 = b2.
  158.    This is the same way that `strcmp' reports the result of comparing
  159.    strings.  */ 
  160.  
  161. int
  162. decimal_compare (b1, b2)
  163.      decimal b1, b2;
  164. {
  165.   int l1, l2;
  166.   char *p1, *p2, *s1, *s2;
  167.   int i;
  168.  
  169.   /* If signs differ, deduce result from the signs */
  170.  
  171.   if (b2->sign && !b1->sign) return 1;
  172.   if (b1->sign && !b2->sign) return -1;
  173.  
  174.   /* If same sign but number of nonfraction digits differs,
  175.      the one with more of them is farther from zero.  */
  176.  
  177.   if (b1->before != b2->before)
  178.     if (b1->sign)
  179.       return (int) (b2->before - b1->before);
  180.     else
  181.       return (int) (b1->before - b2->before);
  182.  
  183.   /* Else compare the numbers digit by digit from high end */
  184.   l1 = LENGTH(b1);
  185.   l2 = LENGTH(b2);  
  186.   s1 = b1->contents;        /* Start of number -- don't back up digit pointer past here */
  187.   s2 = b2->contents;
  188.   p1 = b1->contents + l1;    /* Scanning pointer, for fetching digits.  */
  189.   p2 = b2->contents + l2;
  190.   for (i = MAX(l1, l2); i >= 0; i--)
  191.     {
  192.       int r = ((p1 != s1) ? *--p1 : 0) - ((p2 != s2) ? *--p2 : 0);
  193.       if (r)
  194.     return b1->sign ? -r : r;
  195.     }
  196.   return 0;
  197. }
  198.  
  199. /* Return the number of digits stored in decimal number `b' */
  200.  
  201. int
  202. decimal_length (b)
  203.      decimal b;
  204. {
  205.   return LENGTH(b);
  206. }
  207.  
  208. /* Return the number of fraction digits stored in decimal number `b'.  */
  209.  
  210. int
  211. decimal_after (b)
  212.      decimal b;
  213. {
  214.   return b->after;
  215. }
  216.  
  217. /* Round decimal number `b' to have only `digits' fraction digits.
  218.    Result is rounded to nearest unit in the last remaining digit.
  219.    Return the result, another decimal number.  */
  220.  
  221. decimal
  222. decimal_round_digits (b, digits)
  223.      decimal b;
  224.      int digits;
  225. {
  226.   decimal result;
  227.   int old;
  228.  
  229.   if (b->after <= digits) return decimal_copy (b);
  230.  
  231.   if (digits < 0)
  232.     {
  233.       decimal_error ("request to keep negative number of digits %d", digits);
  234.       return decimal_copy (b);
  235.     }
  236.  
  237.   result = make_decimal (b->before + 1, b->after);
  238.   result->sign = b->sign;
  239.   bcopy (b->contents, result->contents, LENGTH(b));
  240.  
  241.   old = result->after;
  242.  
  243.   /* Add .5 * last place to keep, so that we round rather than truncate */
  244.   /* Note this ignores sign of result, so if result is negative
  245.      it is subtracting */
  246.  
  247.   add_scaled (result, DECIMAL_HALF, 1, old - digits - 1);
  248.  
  249.   /* Flush desired digits, and any trailing zeros exposed by them.  */
  250.  
  251.   flush_trailing_digits (result, old - digits);
  252.  
  253.   /* Flush leading digits -- always is one, unless was a carry into it */
  254.  
  255.   while (result->before > 0
  256.      && result->contents[LENGTH(result) - 1] == 0)
  257.     result->before--;
  258.  
  259.   return result;
  260. }
  261.  
  262. /* Truncate decimal number `b' to have only `digits' fraction digits.
  263.    Any fraction digits in `b' beyond that are dropped and ignored.
  264.    Truncation is toward zero.
  265.    Return the result, another decimal number.  */
  266.  
  267. decimal
  268. decimal_trunc_digits (b, digits)
  269.      decimal b;
  270.      int digits;
  271. {
  272.   decimal result = decimal_copy (b);
  273.   int old = result->after;
  274.  
  275.   if (old <= digits) return result;
  276.  
  277.   if (digits < 0)
  278.     {
  279.       decimal_error ("request to keep negative number of digits %d", digits);
  280.       return result;
  281.     }
  282.  
  283.   flush_trailing_digits (result, old - digits);
  284.  
  285.   return result;
  286. }
  287.  
  288. /* Return the fractional part of decimal number `b':
  289.    that is, `b' - decimal_trunc_digits (`b') */
  290.  
  291. decimal
  292. decimal_fraction (b)
  293.      decimal b;
  294. {
  295.   decimal result = make_decimal (0, b->after);
  296.   bcopy (b->contents, result->contents, b->after);
  297.   return result;
  298. }
  299.  
  300. /* return an integer whose value is that of decimal `b', sans its fraction.  */
  301.  
  302. int
  303. decimal_to_int (b)
  304.      decimal b;
  305. {
  306.   int result = 0;
  307.   int i;
  308.   int end = b->after;
  309.  
  310.   for (i = LENGTH(b) - 1; i >= end; i--)
  311.     {
  312.       result *= RADIX;
  313.       result += b->contents[i];
  314.     }
  315.   return result;
  316. }
  317.  
  318. /* return a decimal whose value is the integer i.  */
  319.  
  320. decimal
  321. decimal_from_int (i)
  322.      int i;
  323. {
  324.   int log, tem;
  325.   decimal result;
  326.  
  327.   for (log = 0, tem = (i > 0 ? i : - i); tem; log++, tem /= RADIX);
  328.  
  329.   result = make_decimal (log, 0);
  330.  
  331.   for (log = 0, tem = (i > 0 ? i : - i); tem; log++, tem /= RADIX)
  332.     result->contents[log] = tem % RADIX;
  333.  
  334.   if (i < 0) result->sign = 1;
  335.   return result;
  336. }
  337.  
  338. /* Return (as an integer) the result of dividing decimal number `b' by
  339.    integer `divisor'. 
  340.    This is used in printing decimal numbers in other radices. */
  341.  
  342. int
  343. decimal_int_rem (b, divisor)
  344.      decimal b;
  345.      int divisor;
  346. {
  347.   int len = LENGTH(b);
  348.   int end = b->after;
  349.   int accum = 0;
  350.   int i;
  351.  
  352.   for (i = len - 1; i >= end; i--)
  353.     {
  354.       accum %= divisor;
  355.       accum *= RADIX;
  356.       accum += b->contents[i];
  357.     }
  358.   return accum % divisor;
  359. }
  360.  
  361. /* Convert digit `digit' to a character and output it by calling
  362.    `charout' with it as arg. */
  363.  
  364. static void
  365. print_digit (digit, charout)
  366.      int digit;
  367.      void (*charout) ();
  368. {
  369.   if (digit < 10)
  370.     charout ('0' + digit);
  371.   else
  372.     charout ('A' + digit - 10);
  373. }
  374.  
  375. /* print decimal number `b' in radix `radix', assuming it is an integer.
  376.    `r' is `radix' expressed as a decimal number. */
  377.  
  378. static
  379. decimal_print_1 (b, r, radix, charout)
  380.      decimal b, r;
  381.      int radix;
  382.      void (*charout) ();
  383. {
  384.   int digit = decimal_int_rem (b, radix);
  385.   decimal rest = decimal_div (b, r, 0);
  386.  
  387.   if (!decimal_zerop (rest))
  388.     decimal_print_1 (rest, r, radix, charout);
  389.  
  390.   print_digit (digit, charout);
  391.  
  392.   free (rest);
  393. }
  394.  
  395. /* User entry: print decimal number `b' in radix `radix' (an integer),
  396.    outputting characters by calling `charout'.  */
  397.  
  398. void
  399. decimal_print (b, charout, radix, precision)
  400.      decimal b;
  401.      void (*charout) ();
  402.      int radix;
  403.      int precision;
  404. {
  405.   if (b->sign) charout ('-');
  406.  
  407.   if (radix == RADIX)
  408.     {
  409.       /* decimal output => just print the digits, inserting a point in
  410.      the proper place.  */ 
  411.       int i;
  412.       int before = b->before;
  413.       int len = before + b->after;
  414.       for (i = 0; i < len; i++)
  415.     {
  416.       if (i == before) charout ('.');
  417.       /* Broken if RADIX /= 10
  418.          charout ('0' + b->contents [len - 1 - i]); */
  419.       print_digit (b->contents [len - 1 - i], charout);
  420.     }
  421.       if (!len)
  422.     charout ('0');
  423.     }
  424.   else
  425.     {
  426.       /* nonstandard radix: must use multiply and divide to determine the
  427.      digits of the number in that radix.  */
  428.  
  429.       int i;
  430.       extern double log10 ();
  431.       /* Compute the number of fraction digits we want to have in the
  432.          new radix.  They should contain the same amount of
  433.          information as the decimal digits we have.  */
  434. /*      int nfrac = (b->after / log10 ((double) radix) + .99); */
  435.       int nfrac = precision;
  436.       decimal r = decimal_from_int (radix);
  437.       decimal intpart = decimal_trunc_digits (b, 0);
  438.  
  439.       /* print integer part */
  440.       decimal_print_1 (intpart, r, radix, charout);
  441.       free (intpart);
  442.  
  443.       /* print fraction part */
  444.       if (nfrac)
  445.     {
  446.           decimal tem1, tem2;
  447.       tem1 = decimal_fraction (b);
  448.       charout ('.');
  449.       /* repeatedly multiply by `radix', print integer part as one digit,
  450.          and flush the integer part.  */
  451.       for (i = 0; i < nfrac; i++)
  452.         {
  453.           if (decimal_zerop(tem1))
  454.         break;
  455.           tem2 = decimal_mul (tem1, r);
  456.           free (tem1);
  457.           print_digit (decimal_to_int (tem2), charout);
  458.           tem1 = decimal_fraction (tem2);
  459.           free (tem2);
  460.         }
  461.       free (tem1);
  462.     }
  463.       free (r);
  464.     }
  465. }
  466.  
  467. static int
  468. decode_digit (digitchar)
  469.      char digitchar;
  470. {
  471.   if ('0' <= digitchar && digitchar <= '9')
  472.     return digitchar - '0';
  473.   if ('a' <= digitchar && digitchar <= 'z')
  474.     return digitchar - 'a' + 10;
  475.   if ('A' <= digitchar && digitchar <= 'Z')
  476.     return digitchar - 'A' + 10;
  477.   return -1;
  478. }
  479.  
  480. /* Parse string `s' into a number using radix `radix'
  481.    and return result as a decimal number.  */
  482.  
  483. decimal
  484. decimal_parse (s, radix, precision)
  485.      char *s;
  486.      int radix;
  487.      int precision;
  488. {
  489.   int i, len, before = -1;
  490.   char *p;
  491.   char c;
  492.   decimal result;
  493.   int negative = 0;
  494.   int excess_digit = 0;
  495.  
  496.   if (*s == '-')
  497.     {
  498.       s++;
  499.       negative = 1;
  500.     }
  501.  
  502.   /* First scan for valid characters.
  503.      Count total num digits, and count num before the decimal point.  */
  504.  
  505.   p = s;
  506.   i = 0;
  507.   while (c = *p++)
  508.     {
  509.       if (c == '.')
  510.         {
  511.       if (before >= 0)
  512.         decimal_error ("two decimal points in %s", s);
  513.           before = i;
  514.     }
  515.       else if (c == '0' && !i && before < 0)
  516.     s++;   /* Discard leading zeros */
  517.       else if (decode_digit (c) >= 0)
  518.     {
  519.       i++;
  520.       if (decode_digit (c) > RADIX)
  521.         excess_digit = 1;
  522.     }
  523.       else
  524.     decimal_error ("invalid number %s", s);
  525.     }
  526.  
  527.   len = i;
  528.   if (before < 0) before = i;
  529.  
  530.   p = s;
  531.  
  532.   /* Now parse those digits */
  533.  
  534.   if (radix != RADIX || excess_digit)
  535.     {
  536.       decimal r = decimal_from_int (radix);
  537.       extern double log10 ();
  538. /*      int digits = (len - before) * log10 ((double) radix) + .99; */
  539.       int digits = precision;
  540.       result = decimal_copy (DECIMAL_ZERO);
  541.  
  542.       /* Parse all the digits into an integer, ignoring decimal point,
  543.      by multiplying by `radix'.  */
  544.  
  545.       while (i > 0 && (c = *p++))
  546.     {
  547.       if (c != '.')
  548.         {
  549.           decimal newdig = decimal_from_int (decode_digit (c));
  550.           decimal prod = decimal_mul (result, r);
  551.           decimal newresult = decimal_add (newdig, prod);
  552.  
  553.           free (newdig);  free (prod);  free (result);
  554.           result = newresult;
  555.           i--;
  556.         }
  557.     }
  558.  
  559.       /* Now put decimal point in right place
  560.      by dividing by `radix' once for each digit
  561.      that really should have followed the decimal point.  */
  562.  
  563.       for (i = before; i < len; i++)
  564.     {
  565.       decimal newresult = decimal_div (result, r, digits);
  566.       free (result);
  567.       result = newresult;
  568.     }
  569.       free (r);
  570.     }
  571.   else
  572.     {
  573.       /* radix is standard - just copy the digits into a decimal number.  */
  574.  
  575.       int tem;
  576.       result = make_decimal (before, len - before);
  577.  
  578.       while (i > 0 && (c = *p++))
  579.     {
  580.       if ((c != '.') &&
  581.           ((tem = decode_digit (c)) >= 0))
  582.         result->contents [--i] = tem;
  583.     }
  584.     }
  585.  
  586.   if (negative) result->sign = 1;
  587.   flush_trailing_digits (result, 0);
  588.   return result;
  589. }
  590.  
  591. /* Add b1 and b2, considering their signs */
  592.  
  593. decimal
  594. decimal_add (b1, b2)
  595.      decimal b1, b2;
  596. {
  597.   decimal v;
  598.  
  599.   if (b1->sign != b2->sign)
  600.     v = decimal_sub1 (b1, b2);
  601.   else
  602.     v = decimal_add1 (b1, b2);
  603.   if (b1->sign && !decimal_zerop (v))
  604.     v->sign = !v->sign;
  605.   return v;
  606. }
  607.  
  608. /* Add b1 and minus b2, considering their signs */
  609.  
  610. decimal
  611. decimal_sub (b1, b2)
  612.      decimal b1, b2;
  613. {
  614.   decimal v;
  615.  
  616.   if (b1->sign != b2->sign)
  617.     v = decimal_add1 (b1, b2);
  618.   else
  619.     v = decimal_sub1 (b1, b2);
  620.   if (b1->sign && !decimal_zerop (v))
  621.     v->sign = !v->sign;
  622.   return v;
  623. }
  624.  
  625. /* Return the negation of b2.  */
  626.  
  627. decimal
  628. decimal_neg (b2)
  629.      decimal b2;
  630. {
  631.   decimal v = decimal_copy (b2);
  632.  
  633.   if (!decimal_zerop (v))
  634.     v->sign = !v->sign;
  635.   return v;
  636. }
  637.  
  638. /* add magnitudes of b1 and b2, ignoring their signs. */
  639.  
  640. static decimal
  641. decimal_add1 (b1, b2)
  642.      decimal b1, b2;
  643. {
  644.   int before = MAX (b1->before, b2->before);
  645.   int after = MAX (b1->after, b2->after);
  646.  
  647.   int len = before+after+1;
  648.   decimal result = make_decimal (before+1, after);
  649.  
  650.   int i;
  651.   char *s1 = b1->contents;
  652.   char *s2 = b2->contents;
  653.   char *p1 = s1 + b1->after - after;
  654.   char *p2 = s2 + b2->after - after;
  655.   char *e1 = s1 + b1->before + b1->after;
  656.   char *e2 = s2 + b2->before + b2->after;
  657.   char *pr = result->contents;
  658.   int accum = 0;
  659.  
  660.   for (i = 0; i < len; i++, p1++, p2++)
  661.     {
  662.       accum /= RADIX;
  663.       if (p1 >= s1 && p1 < e1) accum += *p1;
  664.       if (p2 >= s2 && p2 < e2) accum += *p2;
  665.       *pr++ = accum % RADIX;
  666.     }
  667.   if (!accum)
  668.     (result->before)--;
  669.  
  670.   flush_trailing_digits (result, 0);
  671.  
  672.   return result;
  673. }
  674.  
  675. /* subtract magnitude of b2 from that or b1, returning signed decimal
  676.    number. */ 
  677.  
  678. static decimal
  679. decimal_sub1 (b1, b2)
  680.      decimal b1, b2;
  681. {
  682.   int before = MAX (b1->before, b2->before);
  683.   int after = MAX (b1->after, b2->after);
  684.  
  685.   int len = before+after;
  686.   decimal result = make_decimal (before, after);
  687.  
  688.   int i;
  689.   char *s1 = b1->contents;
  690.   char *s2 = b2->contents;
  691.   char *p1 = s1 + b1->after - after;
  692.   char *p2 = s2 + b2->after - after;
  693.   char *e1 = s1 + b1->before + b1->after;
  694.   char *e2 = s2 + b2->before + b2->after;
  695.   char *pr = result->contents;
  696.   int accum = 0;
  697.  
  698.   for (i = 0; i < len; i++, p1++, p2++)
  699.     {
  700.       if (p1 >= s1 && p1 < e1) accum += *p1;
  701.       if (p2 >= s2 && p2 < e2) accum -= *p2;
  702.       if (accum < 0 && accum % RADIX)
  703.         *pr = RADIX - (- accum) % RADIX;
  704.       else
  705.     *pr = accum % RADIX;
  706.       accum -= *pr++;
  707.       accum /= RADIX;
  708.     }
  709.  
  710.   /* If result is negative, subtract it from RADIX**length
  711.      so that we get the right digits for sign-magnitude
  712.      rather than RADIX-complement */
  713.  
  714.   if (accum)
  715.     {
  716.       result->sign = 1;
  717.       pr = result->contents;
  718.       accum = 0;
  719.       for (i = 0; i < len; i++)
  720.     {
  721.       accum -= *pr;
  722.       if (accum)
  723.         *pr = accum + RADIX;
  724.       else
  725.         *pr = 0;
  726.       accum -= *pr++;
  727.       accum /= RADIX;
  728.     }
  729.     }
  730.  
  731.   /* flush leading nonfraction zero digits */
  732.  
  733.   while (result->before && *--pr == 0)
  734.     (result->before)--;
  735.  
  736.   flush_trailing_digits (result, 0);
  737.  
  738.   return result;
  739. }
  740.  
  741. /* multiply b1 and b2 keeping `digits' fraction digits */
  742.  
  743. decimal
  744. decimal_mul_rounded (b1, b2, digits)
  745.      decimal b1, b2;
  746.      int digits;
  747. {
  748.   decimal tem = decimal_mul (b1, b2);
  749.   decimal result = decimal_round_digits (tem, digits);
  750.   free (tem);
  751.   return result;
  752. }
  753.  
  754. /* multiply b1 and b2 keeping the right number of fraction digits
  755.    for the `dc' program with precision = `digits'.  */
  756.  
  757. decimal
  758. decimal_mul_dc (b1, b2, digits)
  759.      decimal b1, b2;
  760.      int digits;
  761. {
  762.   decimal tem = decimal_mul (b1, b2);
  763.   decimal result
  764.     = decimal_round_digits (tem, MAX (digits, MAX (b1->after, b2->after)));
  765.   free (tem);
  766.   return result;
  767. }
  768.  
  769. /* multiply b1 and b2 as decimal error-free values;
  770.    keep LENGTH(b1) plus LENGTH(b2) significant figures. */
  771.  
  772. decimal
  773. decimal_mul (b1, b2)
  774.      decimal b1, b2;
  775. {
  776.   decimal result = make_decimal (b1->before + b2->before, b1->after + b2->after);
  777.   int i;
  778.   int length2 = LENGTH(b2);
  779.   char *pr;
  780.  
  781.   for (i = 0; i < length2; i++)
  782.     add_scaled (result, b1, b2->contents[i], i);
  783.  
  784.   /* flush leading nonfraction zero digits */
  785.  
  786.   pr = result->contents + LENGTH(result);
  787.   while (result->before && *--pr == 0)
  788.     (result->before)--;
  789.  
  790.   flush_trailing_digits (result, 0);   /* flush trailing zeros */
  791.  
  792.   /* Set sign properly */
  793.  
  794.   if (b1->sign != b2->sign && LENGTH(result))
  795.     result->sign = 1;
  796.  
  797.   return result;
  798. }
  799.  
  800. /* Modify decimal number `into' by adding `from',
  801.    multiplied by `factor' (which should be nonnegative and less than RADIX)
  802.    and shifted left `scale' digits at the least significant end. */
  803.  
  804. static void
  805. add_scaled (into, from, factor, scale)
  806.      decimal into, from;
  807.      int factor, scale;
  808. {
  809.   char *pf = from->contents;
  810.   char *pi = into->contents + scale;
  811.   int lengthf = LENGTH(from);
  812.   int lengthi = LENGTH(into) - scale;
  813.   
  814.   int accum = 0;
  815.   int i;
  816.  
  817.   for (i = 0; i < lengthi; i++)
  818.     {
  819.       accum /= RADIX;
  820.       if (i < lengthf)
  821.         accum += *pf++ * factor;
  822.       accum += *pi;
  823.       *pi++ = accum % RADIX;
  824.     }
  825. }
  826.  
  827.  
  828. /* A utility routine for the divide:  First a one digit multiply.
  829.    NUM (with SIZE digits) is multiplied by DIGIT and the result is
  830.    placed into RESULT.  It is written so that NUM and RESULT can be
  831.    the same pointers.  */
  832.  
  833. static void
  834. _one_mult (num, size, digit, result)
  835.      unsigned char *num;
  836.      int size, digit;
  837.      unsigned char *result;
  838. {
  839.   int carry, value;
  840.   unsigned char *nptr, *rptr;
  841.  
  842.   if (digit == 0)
  843.     memset (result, 0, size);
  844.   else
  845.     {
  846.       if (digit == 1)
  847.         memcpy (result, num, size);
  848.       else
  849.         {
  850.           /* Initialize */
  851.           nptr = (unsigned char *) (num+size-1);
  852.           rptr = (unsigned char *) (result+size-1);
  853.           carry = 0;
  854.  
  855.           while (size-- > 0)
  856.             {
  857.               value = *nptr-- * digit + carry;
  858.               *rptr-- = value % 10;
  859.               carry = value / 10;
  860.             }
  861.  
  862.           if (carry != 0) *rptr = carry;
  863.         }
  864.     }
  865. }
  866.  
  867. /* The full division routine. This computes b1 / b2.  It returns the
  868.    result of the division if ok and DECIMAL_ZERO if not.  The number of
  869.    digits after the decimal point is `digits'.  The algorithm is found
  870.    in Knuth Vol 2. p237. */
  871.  
  872. decimal
  873. decimal_div (b1, b2, digits)
  874.      decimal b1, b2;
  875.      int digits;
  876. {
  877.   decimal n1, n2;
  878.   decimal qval;
  879.   unsigned char *num1, *num2;
  880.   unsigned char *ptr1, *ptr2, *n2ptr, *qptr;
  881.   int  scale1, val;
  882.   unsigned int  len1, len2, scale2, qdigits, extra, count;
  883.   unsigned int  qdig, qguess, borrow, carry;
  884.   unsigned char *mval;
  885.   char zero;
  886.   unsigned int  norm;
  887.  
  888.  
  889.   /* Test for divide by zero. */
  890.   if (decimal_zerop(b2)) {
  891.     decimal_error("divisor is zero", 0);
  892.     return decimal_copy (DECIMAL_ZERO);
  893.   }
  894.  
  895.   /* Test for divide by 1.  If it is we must truncate. */
  896.   if (b2->after == 0)
  897.     {
  898.       if (b2->before == 1 && *b2->contents == 1)
  899.         {
  900.       int size_to_copy, copy_to_offset, copy_from_offset;
  901.  
  902.           qval = make_decimal(b1->before, digits);
  903.           qval->sign = (b1->sign == b2->sign ? 0 : 1);
  904.           memset (&qval->contents[b1->before], 0, digits);
  905.  
  906.       size_to_copy = b1->before + MIN(b1->after, digits);
  907.       copy_to_offset = qval->before + qval->after - size_to_copy;
  908.           memcpy (&qval->contents[copy_to_offset], b1->contents + MAX(b1->after-digits, 0), b1->before + MIN(b1->after,digits));
  909.       flush_trailing_digits(qval, 0);
  910.       return qval;
  911.         }
  912.     }
  913.  
  914.   /* reverse b1 and b2 */
  915.   n1 = decimal_copy(b1);
  916.   n2 = decimal_copy(b2);
  917.   {
  918.     int i, len;
  919.  
  920.     len = n1->before + n1->after;
  921.     for (i = 0; i < len; i++)
  922.         n1->contents[len - i - 1] = b1->contents[i];
  923.  
  924.     len = n2->before + n2->after;
  925.     for (i = 0; i < len; i++)
  926.         n2->contents[len - i - 1] = b2->contents[i];
  927.   }
  928.   
  929.  
  930.   /* Set up the divide.  Move the decimal point on n1 by n2's scale.
  931.      Remember, zeros on the end of num2 are wasted effort for dividing. */
  932.   scale2 = n2->after;
  933.   n2ptr = (unsigned char *) n2->contents + n2->before + scale2 - 1;
  934.   while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
  935.  
  936.   len1 = n1->before + scale2;
  937.   scale1 = n1->after - scale2;
  938.   if (scale1 < digits)
  939.     extra = digits - scale1;
  940.   else
  941.     extra = 0;
  942.   num1 = (unsigned char *) malloc (n1->before+n1->after+extra+2);
  943.   if (!num1) {
  944.     decimal_error("Fatal error: Out of memory for malloc.\n", 0);
  945.     exit(1);
  946.   }
  947.   memset (num1, 0, n1->before+n1->after+extra+2);
  948.   memcpy (num1+1, n1->contents, n1->before+n1->after);
  949.  
  950.   len2 = n2->before + scale2;
  951.   num2 = (unsigned char *) malloc (len2+1);
  952.   if (!num2) {
  953.     decimal_error("Fatal error: Out of memory for malloc.\n", 0);
  954.     exit(1);
  955.   }
  956.   memcpy (num2, n2->contents, len2);
  957.   *(num2+len2) = 0;
  958.   n2ptr = num2;
  959.   while (*n2ptr == 0)
  960.     {
  961.       n2ptr++;
  962.       len2--;
  963.     }
  964.  
  965.   /* Calculate the number of quotient digits. */
  966.   if (len2 > len1+digits)
  967.     {
  968.       qdigits = digits+1;
  969.       zero = 1;
  970.     }
  971.   else
  972.     {
  973.       zero = 0;
  974.       if (len2>len1)
  975.         qdigits = digits+1;      /* One for the zero integer part. */
  976.       else
  977.         qdigits = len1-len2+digits+1;
  978.     }
  979.  
  980.   /* Allocate and zero the storage for the quotient. */
  981.   qval = make_decimal (qdigits-digits,digits);
  982.   memset (qval->contents, 0, qdigits);
  983.  
  984.   /* Allocate storage for the temporary storage mval. */
  985.   mval = (unsigned char *) malloc (len2+1);
  986.   if (!mval) {
  987.     decimal_error("Fatal error: Out of memory for malloc.\n", 0);
  988.     exit(1);
  989.   }
  990.  
  991.   /* Now for the full divide algorithm. */
  992.   if (!zero)
  993.     {
  994.       /* Normalize */
  995.       norm =  10 / ((int)*n2ptr + 1);
  996.       if (norm != 1)
  997.         {
  998.           _one_mult (num1, len1+scale1+extra+1, norm, num1);
  999.           _one_mult (n2ptr, len2, norm, n2ptr);
  1000.         }
  1001.  
  1002.       /* Initialize divide loop. */
  1003.       qdig = 0;
  1004.       if (len2 > len1)
  1005.         qptr = (unsigned char *) qval->contents+len2-len1;
  1006.       else
  1007.         qptr = (unsigned char *) qval->contents;
  1008.  
  1009.       /* Loop */
  1010.       while (qdig <= len1+digits-len2)
  1011.         {
  1012.           /* Calculate the quotient digit guess. */
  1013.           if (*n2ptr == num1[qdig])
  1014.             qguess = 9;
  1015.           else
  1016.             qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
  1017.  
  1018.           /* Test qguess. */
  1019.           if (n2ptr[1]*qguess >
  1020.               (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
  1021.                + num1[qdig+2])
  1022.             {
  1023.               qguess--;
  1024.               /* And again. */
  1025.               if (n2ptr[1]*qguess >
  1026.                   (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
  1027.                   + num1[qdig+2])
  1028.                 qguess--;
  1029.             }
  1030.  
  1031.           /* Multiply and subtract. */
  1032.           borrow = 0;
  1033.           if (qguess != 0)
  1034.             {
  1035.               *mval = 0;
  1036.               _one_mult (n2ptr, len2, qguess, mval+1);
  1037.               ptr1 = (unsigned char *) num1+qdig+len2;
  1038.               ptr2 = (unsigned char *) mval+len2;
  1039.               for (count = 0; count < len2+1; count++)
  1040.                 {
  1041.                   val = (int) *ptr1 - (int) *ptr2-- - borrow;
  1042.                   if (val < 0)
  1043.                     {
  1044.                       val += 10;
  1045.                       borrow = 1;
  1046.                     }
  1047.                   else
  1048.                     borrow = 0;
  1049.                   *ptr1-- = val;
  1050.                 }
  1051.             }
  1052.  
  1053.           /* Test for negative result. */
  1054.           if (borrow == 1)
  1055.             {
  1056.               qguess--;
  1057.               ptr1 = (unsigned char *) num1+qdig+len2;
  1058.               ptr2 = (unsigned char *) n2ptr+len2-1;
  1059.               carry = 0;
  1060.               for (count = 0; count < len2; count++)
  1061.                 {
  1062.                   val = (int) *ptr1 + (int) *ptr2-- + carry;
  1063.                   if (val > 9)
  1064.                     {
  1065.                       val -= 10;
  1066.                       carry = 1;
  1067.                     }
  1068.                   else
  1069.                     carry = 0;
  1070.                   *ptr1-- = val;
  1071.                 }
  1072.               if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
  1073.             }
  1074.  
  1075.           /* We now know the quotient digit. */
  1076.           *qptr++ =  qguess;
  1077.           qdig++;
  1078.         }
  1079.     }
  1080.  
  1081.   /* must reverse qval */
  1082.   {
  1083.     int i, len;
  1084.     decimal temp = decimal_copy(qval);
  1085.  
  1086.     len = qval->before + qval->after;
  1087.     for (i = 0; i < len; i++)
  1088.         qval->contents[len - i - 1] = temp->contents[i];
  1089.     free (temp);
  1090.   }
  1091.  
  1092.   /* Clean up and return the number... */
  1093.  
  1094.   /* first flush trailing zero fraction digits */
  1095.   flush_trailing_digits (qval, 0);
  1096.  
  1097.   /* then flush leading nonfraction zero digits */
  1098.   {
  1099.     char *pr = qval->contents + qval->before + qval->after;
  1100.     while (qval->before && *--pr == 0)
  1101.       (qval->before)--;
  1102.   }
  1103.  
  1104.   /* finally fix up the sign */
  1105.   qval->sign = (n1->sign != n2->sign);
  1106.   if (decimal_zerop (qval)) qval->sign = 0;
  1107.  
  1108.   /* Clean up temporary storage. */
  1109.   free (mval);
  1110.   free (num1);
  1111.   free (num2);
  1112.   free (n1);
  1113.   free (n2);
  1114.  
  1115.   return qval;
  1116. }
  1117.  
  1118. /* The remainder from `b1' / `b2' is equivalent to:
  1119.    `b1' - (`b1' / `b2') * `b2'.
  1120. */
  1121.  
  1122. decimal
  1123. decimal_rem(b1, b2, digits)
  1124.      decimal b1, b2;
  1125.      int digits;
  1126. {
  1127.   decimal t1, t2, t3;
  1128.   int rscale;
  1129.  
  1130.   /* return error if division by zero attempted */
  1131.   if (decimal_zerop(b2)) {
  1132.     decimal_error("divisor is zero", 0);
  1133.     return decimal_copy(DECIMAL_ZERO);
  1134.   }
  1135.  
  1136.   /* calculate final scale...max of b1->after and b2->after+digits */
  1137.   rscale = b1->after > b2->after+digits ? b1->after : b2->after+digits;
  1138.  
  1139.   /* calculate answer */
  1140.   t1 = decimal_div(b1, b2, digits);
  1141.   t2 = decimal_mul(t1, b2, rscale);
  1142.   t3 = decimal_sub(b1, t2);
  1143.   free(t1);
  1144.   free(t2);
  1145.  
  1146.   return t3;
  1147. }
  1148.  
  1149.  
  1150. /* returns negative number if we chose factor too large */
  1151.  
  1152. static int
  1153. subtract_scaled (into, from, factor, scale)
  1154.      decimal into, from;
  1155.      int factor, scale;
  1156. {
  1157.   char *pf = from->contents;
  1158.   char *pi = into->contents + scale;
  1159.   int lengthf = LENGTH(from);
  1160.   int lengthi = LENGTH(into) - scale;
  1161.   int accum = 0;
  1162.   int i;
  1163.  
  1164.   for (i = 0; i < lengthi && i <= lengthf; i++)
  1165.     {
  1166.       if (i < lengthf)
  1167.         accum -= *pf++ * factor;
  1168.       accum += *pi;
  1169.       if (accum < 0 && accum % RADIX)
  1170.         *pi = RADIX - (- accum) % RADIX;
  1171.       else
  1172.     *pi = accum % RADIX;
  1173.       accum -= *pi++;
  1174.       accum /= RADIX;
  1175.     }
  1176.   return accum;
  1177. }
  1178.  
  1179. /* Return the square root of decimal number D, using Newton's method.
  1180.    Number of fraction digits returned is max of FRAC_DIGITS
  1181.    and D's number of fraction digits.  */
  1182.  
  1183. decimal
  1184. decimal_sqrt (d, frac_digits)
  1185.      decimal d;
  1186.      int frac_digits;
  1187. {
  1188.   decimal guess;
  1189.   int notdone = 1;
  1190.  
  1191.   if (decimal_zerop (d)) return d;
  1192.   if (d->sign)
  1193.     {
  1194.       decimal_error ("square root argument negative", 0);
  1195.       return decimal_copy (DECIMAL_ZERO);
  1196.     }
  1197.  
  1198.   frac_digits = MAX (frac_digits, d->after);
  1199.  
  1200.   /* Compute an initial guess by taking the square root 
  1201.      of a nearby power of RADIX.  */
  1202.  
  1203.   if (d->before)
  1204.     {
  1205.       guess = make_decimal ((d->before + 1) / 2, 0);
  1206.       guess->contents[guess->before - 1] = 1;
  1207.     }
  1208.   else
  1209.     {
  1210.       /* Arg is less than 1; compute nearest power of RADIX */
  1211.       char *p = d->contents + LENGTH(d);
  1212.       char *sp = p;
  1213.  
  1214.       while (!*--p);    /* Find most significant nonzero digit */
  1215.       if (sp - p == 1)
  1216.     {
  1217.       /* Arg is bigger than 1/RADIX; use 1 as a guess */
  1218.       guess = decimal_copy (DECIMAL_ONE);
  1219.     }
  1220.       else
  1221.     {
  1222.       guess = make_decimal (0, (sp - p) / 2);
  1223.       guess->contents[0] = 1;
  1224.     }
  1225.     }
  1226.  
  1227.   /* Iterate doing guess = (guess + d/guess) / 2  */
  1228.  
  1229.   while (notdone)
  1230.     {
  1231.       decimal tem1 = decimal_div (d, guess, frac_digits + 1);
  1232.       decimal tem2 = decimal_add (guess, tem1);
  1233.       decimal tem3 = decimal_mul_rounded (tem2, DECIMAL_HALF, frac_digits);
  1234.       notdone = decimal_compare (guess, tem3);
  1235.       free (tem1);
  1236.       free (tem2);
  1237.       free (guess);
  1238.       guess = tem3;
  1239.       if (decimal_zerop (guess)) return guess;  /* Avoid divide-by-zero */
  1240.     }
  1241.  
  1242.   return guess;
  1243. }
  1244.  
  1245. /* Raise decimal number `base' to power of integer part of decimal
  1246.    number `expt'.
  1247.    This function depends on using radix 10.
  1248.    It is too hard to write it to work for any value of RADIX,
  1249.    so instead it is simply not available if RADIX is not ten.  */
  1250.  
  1251. #if !(RADIX - 10)
  1252.  
  1253. decimal
  1254. decimal_expt (base, expt, frac_digits)
  1255.      decimal base, expt;
  1256.      int frac_digits;
  1257. {
  1258.   decimal accum = decimal_copy (DECIMAL_ONE);
  1259.   decimal basis1 = base;
  1260.   int digits = expt->before;
  1261.   int dig = 0;                /* Expt digit being processed */
  1262.  
  1263.   if (expt->sign)
  1264.   /* If negative power, take reciprocal first thing
  1265.      so that fraction digit truncation won't destroy
  1266.      what will ultimately be nonfraction digits.  */
  1267.     basis1 = decimal_div (DECIMAL_ONE, base, frac_digits);
  1268.   while (dig < digits)
  1269.     {
  1270.       decimal basis2, basis4, basis8, basis10;
  1271.       int thisdigit = expt->contents[expt->after + dig];
  1272.  
  1273.       /* Compute factors to multiply in for each bit of this digit */
  1274.  
  1275.       basis2 = decimal_mul_rounded (basis1, basis1, frac_digits);
  1276.       basis4 = decimal_mul_rounded (basis2, basis2, frac_digits);
  1277.       basis8 = decimal_mul_rounded (basis4, basis4, frac_digits);
  1278.  
  1279.       /* Now accumulate the factors this digit value selects */
  1280.  
  1281.       if (thisdigit & 1)
  1282.     {
  1283.       decimal accum1 = decimal_mul_rounded (accum, basis1, frac_digits);
  1284.       free (accum);
  1285.       accum = accum1;
  1286.     }
  1287.  
  1288.       if (thisdigit & 2)
  1289.     {
  1290.       decimal accum1 = decimal_mul_rounded (accum, basis2, frac_digits);
  1291.       free (accum);
  1292.       accum = accum1;
  1293.     }
  1294.  
  1295.       if (thisdigit & 4)
  1296.     {
  1297.       decimal accum1 = decimal_mul_rounded (accum, basis4, frac_digits);
  1298.       free (accum);
  1299.       accum = accum1;
  1300.     }
  1301.  
  1302.       if (thisdigit & 8)
  1303.     {
  1304.       decimal accum1 = decimal_mul_rounded (accum, basis8, frac_digits);
  1305.       free (accum);
  1306.       accum = accum1;
  1307.     }
  1308.  
  1309.       /* If there are further digits, compute the basis1 for the next digit */
  1310.  
  1311.       if (++dig < digits)
  1312.     basis10 = decimal_mul_rounded (basis2, basis8, frac_digits);
  1313.  
  1314.       /* Free intermediate results */
  1315.  
  1316.       if (basis1 != base) free (basis1);
  1317.       free (basis2);
  1318.       free (basis4);
  1319.       free (basis8);
  1320.       basis1 = basis10;
  1321.     }
  1322.   return accum;
  1323. }
  1324. #endif
  1325.  
  1326. #ifdef TEST
  1327.  
  1328. fputchar (c)
  1329.      char c;
  1330. {
  1331.   putchar (c);
  1332. }
  1333.  
  1334. /* Top level that can be used to test the arithmetic functions */
  1335.  
  1336. main ()
  1337. {
  1338.   char s1[40], s2[40];
  1339.   decimal b1, b2, b3;
  1340.   char c;
  1341.  
  1342.   while (1)
  1343.     {
  1344.       scanf ("%s %c %s", s1, &c, s2);
  1345.       b1 = decimal_parse (s1, RADIX);
  1346.       b2 = decimal_parse (s2, RADIX);
  1347.       switch (c)
  1348.     {
  1349.     default:
  1350.       c = '+';
  1351.     case '+':
  1352.       b3 = decimal_add (b1, b2);
  1353.       break;
  1354.     case '*':
  1355.       b3 = decimal_mul (b1, b2);
  1356.       break;
  1357.         case '/':
  1358.       b3 = decimal_div (b1, b2, 3);
  1359.       break;
  1360.     case '%':
  1361.       b3 = decimal_rem (b1, b2, 3);
  1362.       break;
  1363.         case 'p':
  1364.       decimal_print (b1, fputchar, RADIX);
  1365.       printf (" printed in base %d is ", decimal_to_int (b2));
  1366.       decimal_print (b1, fputchar, decimal_to_int (b2));
  1367.       printf ("\n");
  1368.       continue;
  1369.     case 'r':
  1370.       printf ("%s read in base %d is ", s1, decimal_to_int (b2));
  1371.       decimal_print (decimal_parse (s1, decimal_to_int (b2)), fputchar, RADIX);
  1372.       printf ("\n");
  1373.       continue;
  1374.     }
  1375.       decimal_print (b1, fputchar, RADIX);
  1376.       printf (" %c ", c);
  1377.       decimal_print (b2, fputchar, RADIX);
  1378.       printf (" = ");
  1379.       decimal_print (b3, fputchar, RADIX);
  1380.       printf ("\n");
  1381.     }
  1382. }
  1383.  
  1384. decimal_error (s1, s2)
  1385.      char *s1, *s2;
  1386. {
  1387.   printf ("\n");
  1388.   printf (s1, s2);
  1389.   printf ("\n");
  1390. }
  1391.  
  1392. static void
  1393. pbi (b)
  1394.       int b;
  1395. {
  1396.   decimal_print ((decimal) b, fputchar, RADIX);
  1397. }
  1398.  
  1399. static void
  1400. pb (b)
  1401.       decimal b;
  1402. {
  1403.   decimal_print (b, fputchar, RADIX);
  1404. }
  1405.  
  1406. #endif
  1407.